# JAGS/BUGS Code ----------------------------------------------------------

# n = n.legislators
# m = n.items

# Original model as used in paper
jags_model <- '
model {
	for (i in 1:n) {
		for (j in 1:m) {
			y[i,j] ~ dcat(p[i,j,1:3])
			p[i,j,1] <- Q[i,j,1]
			p[i,j,2] <- Q[i,j,2] - Q[i,j,1]
			p[i,j,3] <- 1 - Q[i,j,2] 
			for(k in 1:2) {
				probit(Q[i,j,k]) <- cut[j, k] - beta[j]*theta[i] 
			}
		}
	}
	
	# Priors
	for(j in 1:m) { beta[j]  ~ dnorm(0, 0.25) }
	for(i in 1:n) { theta[i] ~ dnorm(0, 1) }
	
	for(j in 1:m) { 
		cut[j, 1] ~ dnorm(0, .15) 
		delta[j] ~ dexp(.5)
		cut[j, 2] <- cut[j, 1] + delta[j]
	}
}
'

# Indifference model with actor-specific cutting points
jags_model_actor_cut <- '
model {
  for (i in 1:n) {
    cut[i,1] <- cut[i,2] * (-1)
		for (j in 1:m) {
			y[i,j] ~ dcat(p[i,j,1:3])
			p[i,j,1] <- Q[i,j,1]
			p[i,j,2] <- Q[i,j,2] - Q[i,j,1]
			p[i,j,3] <- 1 - Q[i,j,2] 
			for(k in 1:2) {
				probit(Q[i,j,k]) <- cut[i, k] - beta[j]*theta[i] + alpha[j]
			}
		}
	}
	
	# Priors
  for(j in 1:m) { alpha[j]  ~ dnorm(0, 0.25) }	
  for(j in 1:m) { beta[j]  ~ dnorm(0, 0.25) }
	
  for(i in 1:(n)) { 
  theta[i] ~ dnorm(0, 1) 
  }


	for(j in 1:n) { 
#		cut[j, 1] ~ dnorm(0, .15) 
#   delta[j] ~ dexp(.5)
		cut[j, 2] ~ dunif(0,5)
	}
}
'

# Original model with user-specified priors
jags_model_priors <- '
model {
  for (i in 1:n) {
    for (j in 1:m) {
      y[i,j] ~ dcat(p[i,j,1:3])
      p[i,j,1] <- Q[i,j,1]
      p[i,j,2] <- Q[i,j,2] - Q[i,j,1]
      p[i,j,3] <- 1 - Q[i,j,2] 
      for(k in 1:2) {
        probit(Q[i,j,k]) <- cut[j, k] - beta[j]*theta[i]
      }
    }
  }
  
  # Priors
  for(j in 1:m) { beta[j]  ~ dnorm(0, 0.25) }
  for(i in 1:n) { theta[i] ~ dnorm(t0[i], T0[i]) } 
  
  for(j in 1:m) { 
    cut[j, 1] ~ dnorm(0, .15) 
    delta[j] ~ dexp(.5)
    cut[j, 2] <- cut[j, 1] + delta[j]
  }
}
'

writeLines(jags_model, con="jags_model.bug")
writeLines(jags_model_actor_cut, con="jags_model_actor_cut.bug")
writeLines(jags_model_priors, con="jags_model_priors.bug")


# Install package rjags if not available ----------------------------------
if(!("rjags" %in% library()$results[,1])) {
  cat("Package rjags is not available. Make sure JAGS is installed on your system.\n")
  cat("Download via: http://sourceforge.net/projects/mcmc-jags/files/JAGS/ \n")
  cat("Package rjags is not installed. Installing. \n")
  install.packages("rjags")
}

# Confrontational ---------------------------------------------------------

confrontational = function(dat,
                           priors=NULL,
                           n.burnin=20000,
                           n.samples=100000,
                           n.thin=ifelse(n.samples>500,round(n.samples/10000),1),
                           model.file=NULL,
                           seed=runif(1,0,1e10)) {
  require(rjags, quiet=TRUE)
  
  # Check data
  if (!(is.data.frame(dat)|is.matrix(dat))) stop("Data appears not to be a data frame or matrix")
  values <- unique(unlist(as.data.frame(dat)))
  if(length(values)>3) stop("Data contains more than three unique values.")
  if(!all(values %in% c(1,2,3,NA) & length(values) <= 3)) {
    cat("Data contains values other than 1, 2 and 3. Data will be assumed to be in ordinal form and will be recoded.\n")
    value.sort <- sort(values)
    value.cells <- list()
    for(i in 1:3) {
      value.cells[[i]] <- which(dat==value.sort[i], arr.ind=TRUE)
    }
    for(i in 1:3) {
      dat[value.cells[[i]]] <- i
    }
    cat("Header of recoded dataframe:\n")
    print(head(dat))
    values <- unique(unlist(as.data.frame(dat)))
  }
  if(!all(values %in% c(1,2,3,NA))) stop("Data contains values other than 1,2, or 3. Please recode.")
  
  # Run model in JAGS
  if(is.null(priors)) {
    # Without priors
    if(is.null(model.file)) model.file = "jags_model.bug"
    
    # Set inits to sum model
    thetaInit <- as.numeric(scale(rowSums(dat,na.rm=TRUE)))
    
    model = jags.model(model.file, data=list(y=dat,n=dim(dat)[1],m=dim(dat)[2]),
                       inits=list(.RNG.name="base::Mersenne-Twister", 
                                  .RNG.seed=seed,
                                  theta=thetaInit))
  }  else 
    {
    # With priors
    # Standardize priors, because parties will be estimated on standardized scale
    priorMean <- scale(priors$mean)
    priorPrec <- (1 / scale(sapply((priors$se^2), function(q) ifelse(q==0,.1,q)),  scale=attr(priorMean, "scaled:scale"), center=FALSE))
    
    # Set inits to sum model
    thetaInit <- as.numeric(scale(rowSums(dat,na.rm=TRUE)))
    
    forJags = list(y=dat,n=dim(dat)[1],
                   m=dim(dat)[2],
                   t0=as.numeric(priorMean),
                   T0=as.numeric(priorPrec))
    if(is.null(model.file)) model.file = "jags_model_priors.bug"
    model = jags.model(model.file, forJags,
                       inits=list(.RNG.name="base::Mersenne-Twister", 
                                  .RNG.seed=seed,
                                  theta=thetaInit))
  }
  update(model, n.burnin)
  samples = coda.samples(model, c("theta", "beta", "cut"), n.samples, thin=n.thin)
  
  # Obtain values for theta (party positions), beta (item positions) and tau(cut values)
  theta = samples[[1]][,grep("^theta",dimnames(samples[[1]])[[2]])]
  beta = samples[[1]][,grep("^beta",dimnames(samples[[1]])[[2]])]
  tau = samples[[1]][,grep("^cut",dimnames(samples[[1]])[[2]])]
  
  # Normalize positions
  thetamean = apply(theta, 1, function(x) mean(x))
  thetasd = apply(theta, 1, function(x) sd(x))
  
  thetaStar = apply(theta, 1, function(x) (x-mean(x))/sd(x))
  
  betaStar = beta  
  for(iter in 1:dim(beta)[1]){
    betaStar[iter,] <- beta[iter,]*thetasd[iter]
  }
  
  tauStar = tau
  for(iter in 1:dim(tau)[1]){
    tauStar[iter,] <- tau[iter,]-beta[iter,]*thetamean[iter]
  }
  
  rownames(thetaStar) = rownames(dat)
  colnames(betaStar) = colnames(dat)
  
  # Flip the scale if correlation with summative scores < 0
  # (for ease of interpretation)
  flip <- mean(cor(thetaStar, thetaInit))< 0
  if(flip) {
    thetaStar <- thetaStar * -1
    betaStar <- betaStar * -1
    tauStar <- tauStar * -1
  }
  
  # Create vector with results
  thetaStar.Results = t(apply(thetaStar, 1, function(x)c(mean(x),quantile(x,c(0.025,.975)))))
  betaStar.Results = t(apply(betaStar, 2, function(x)c(mean(x),quantile(x,c(0.025,.975)))))
  tauStar.Results = t(apply(tauStar, 2, function(x)c(mean(x),quantile(x,c(0.025,.975)))))
  
  rownames(thetaStar.Results) = rownames(dat)
  rownames(betaStar.Results) = colnames(dat)
  cn <- c("Mean", "2.5%", "97.5%")
  colnames(thetaStar.Results) <- colnames(betaStar.Results) <- colnames(tauStar.Results) <- cn
  
  
  returnObj <- list(actors = as.mcmc(t(thetaStar)), 
                    items = as.mcmc(betaStar),
                    cutoffs = as.mcmc(tauStar),
                    actors_summary=thetaStar.Results,
                    items_summary=betaStar.Results, 
                    cutoffs_summary=tauStar.Results)
  class(returnObj) <- "confrontational"
  return(returnObj)
}


# summary.confrontational -------------------------------------------------

summary.confrontational <- function(x) {
  cat("Confrontational IRT Model\n")
  cat("-------------------------\n")
  cat("Number of actors: ", nrow(x$actors_summary), "\n")
  cat("Number of items: ", nrow(x$items_summary), "\n")
  cat("\n")
  cat("Actor positions\n")
  print(x$actors_summary)
  
  cat("\n")
  cat("Item positions\n")
  print(x$items_summary)
  return(list(actors_summary=x$actors_summary, 
              items_summary=x$items_summary))
}



# plot.confrontational ----------------------------------------------------

plot.confrontational = function(x,
                                type="actors",
                                labelcex=0.8,
                                main=NULL) {
  
  if(type=="actors") o <- x$actors_summary
  if(type=="items") o <- x$items_summary
  if(type=="cutoffs") o <- x$cutoffs_summary
  
  n = dim(o)[1]
  par(cex.axis=.75,tcl=-0.2,mgp=c(1,.35,0))
  
  plot( x=c(floor(min(o[,2])*1.05),ceiling(max(o[,3])*1.05)),
        y=c(0.25,n+.75),
        yaxs="i",
        xlab="",
        ylab="",
        type="n",
        yaxt="n",
        axes=FALSE, main=main)

  indx = order(o[,1])
  segments(x0=o[indx,2],
           x1=o[indx,3],
           y0=1:n,y1=1:n)
  points(x=o[indx,1],
         y=1:n,
         cex=0.75,
         col="red",
         pch=16,
         xpd=NA)
  ## label
  for(i in 1:n) {
    if(o[indx[i],1]<0) {
      text(y=i,x=o[indx[i],3],
           pos=4,cex=labelcex,labels=dimnames(o)[[1]][indx[i]]) }
    else {
      text(y=i,x=o[indx[i],2],
           pos=2,cex=labelcex,dimnames(o)[[1]][indx[i]])     
    }
  }
  axis(1)
  axis(3)
  abline(v=0, lty=2, col="grey")
}



# biplot.confrontational --------------------------------------------------

biplot.confrontational = function(x, y,
                                type="actors",
                                labelcex=0.8,
                                adj=c(-.15,-.15),
                                main=NULL,
                                labelx="X",
                                labely="Y") {
    
  if(type=="actors") {
    d1 <- x$actors_summary
    d2 <- y$actors_summary
  } 
  if(type=="items") {
    d1 <- x$items_summary
    d2 <- y$items_summary
  }
  if(type=="cutoffs") {
    d1 <- x$cutoffs_summary
    d2 <- y$cutoffs_summary
  }
  
  colnames(d1) <- paste("X", c("m","l","h"), sep=".")
  colnames(d2) <- paste("Y", c("m","l","h"), sep=".")
  
  d <- merge(d1, d2, by=0, all=TRUE)
  colnames(d)[1] <- "actor"
  
  par(cex.axis=.75,tcl=-0.2,mgp=c(1,.35,0))
  
  plot( x=c(floor(min(d[,3])*1.05),ceiling(max(d[,4])*1.05)),
        y=c(floor(min(d[,6])*1.05),ceiling(max(d[,7])*1.05)),
        yaxs="i",
        xlab=labelx,
        ylab=labely,
        type="n",
        yaxt="n",
        axes=FALSE, main=main)
  
  indx = order(d[,2])
  segments(x0=d[indx,3],
           x1=d[indx,4],
           y0=d[indx,5],
           y1=d[indx,5], col="grey50")
  segments(x0=d[indx,2],
           x1=d[indx,2],
           y0=d[indx,6],
           y1=d[indx,7], col="grey50")
  
  colours <- c("red", "blue", "darkgreen", "magenta", "brown")
  
  points(x=d[indx,2],
         y=d[indx,5],
         cex=0.75,
         col=colours,
         pch=16,
         xpd=NA)
  require(plotrix)
  
  thigmophobe.labels(x=d[indx,2],
       y=d[indx,5]+.1,
       labels=d[indx,1],
       cex=labelcex,
       col=colours,
       pch=16,
       xpd=NA)
  
#   text(x=d[indx,2],
#          y=d[indx,5],
#          labels=d[indx,1],
#          cex=labelcex,
#          col=colours,
#          adj=adj,
#          pch=16,
#          xpd=NA)
  axis(2)
  axis(1)
  box()
}



# traceplot.confrontational -----------------------------------------------

traceplot.default <- function(x) {}

traceplot <- function(x, ...)
  UseMethod("traceplot")

traceplot.confrontational = function(x,
                                     type="actors") {
  
  if(type=="actors") dat = x$actors
  if(type=="items") dat = x$items
  if(type=="cutoffs") dat = x$cutoffs
  plot(dat)
}